function results = EstimateSAR(Data,AllW, VarNames)

%% Create table with all data elements
idpair = VarNames.idpair;
wscore = VarNames.wscore;
idvar  = VarNames.idvar ;
tvar   = VarNames.tvar  ;
xvars  = VarNames.xvars ;
yvar   = VarNames.yvar  ;

% Define a cutoff for W
cutoff = 600;

% number of x variables
p = size(xvars,2);


%==================================================================
%           Read/Create the W matrices, and necessary data
%==================================================================

% merge (join) tables to ensure W matches data
ids = unique(AllW(:,[idpair(1),tvar]), 'rows');
ids.Properties.VariableNames{idpair(1)} = idvar;
ids = sortrows(ids,{tvar,idvar});

ModData = sortrows(innerjoin(Data, ids),{tvar,idvar});

% eliminate missing values
TF = (sum(isnan(ModData{:,[yvar, xvars]}), 2) == 0);
ModData    = ModData(TF,:);

% Initialize cell arrays to hold the matrices
begYear = min(ModData{:,tvar});
endYear = max(ModData{:,tvar});

W = cell(1,endYear - begYear + 1);
IDS = cell(1,endYear - begYear + 1);
X = cell(1,endYear - begYear + 1);
Y = cell(1,endYear - begYear + 1);

ID = ModData{:,{idvar,tvar}};

%% process each year of the data
noiseparm = 0.0001;
for i = 1:(endYear - begYear + 1);

    year = begYear + i - 1;

    yeardata = ModData(ModData{:,tvar}==year,:);
    yearW = AllW(AllW{:,tvar}==year,[idpair,tvar,wscore]);

    % keep IDs and generate a new index for IDs
    ids = sortrows(yeardata(:,idvar));
    ids = [ids array2table((1:size(ids,1))')];
    ids.Properties.VariableNames = {idvar,'altid'};
    IDS{i} = ids;
    % Create annual W matrix
    yearW1 = innerjoin(yearW,ids,'LeftKey',1,'RightKey',1);
    yearW1.Properties.VariableNames{'altid'} = 'altid1';

    yearW2 = innerjoin(yearW1,ids,'LeftKey',2,'RightKey',1);
    yearW2.Properties.VariableNames{'altid'} = 'altid2';

    % make sure the matrix is symmetric (some firms might drop)
    iota = ones(size(ids,1),1);
    tmp = array2table([ids{:,'gvkey'}, ids{:,'gvkey'}, year*iota, -noiseparm*iota, ids{:,'altid'}, ids{:,'altid'}]);
    tmp.Properties.VariableNames = yearW2.Properties.VariableNames;
    yearW2 = [yearW2; tmp];

    SCORE = yearW2{:,wscore}  + noiseparm;
    FIRM1 = yearW2{:,'altid1'};
    FIRM2 = yearW2{:,'altid2'};

    W{i} = sparse(FIRM1, FIRM2, SCORE);

    % optionally, restrict W based on cutoff and rescale W
%     n = size(W{i},1);
%     iota = ones(n,1);
%
%     sortedW = sort(W{i},2, 'DESCEND');
% 	W{i}(W{i} < repmat(sortedW(:,cutoff),1,n)) = 0;
    % W{i}(W{i} < prctile(SCORE,cutoff)) = 0;

    %row normalize W
    W{i} = normw(W{i});

    W{i} = sparse(W{i});

    % get X and Y
    Y{i} = studentize(yeardata{:,yvar});
    % Y{i} = yeardata{:,yvar};
    X{i} = studentize(yeardata{:,xvars});
    % X{i} = yeardata{:,xvars};

end;


%==================================================================
%           Stack data to get the proper format
%==================================================================

% process each year of the data
for i = 1:(endYear - begYear + 1);
    year = begYear + i - 1;


    nvec(i,1) = length(Y{i});

	if i == 1
		Wy = W{i} * Y{i};

		yvec = Y{i};

		x = X{i};
	else

		Wy = [ Wy; W{i} * Y{i}];

		yvec = [yvec; Y{i}];

		x = [x; X{i}];

	end
end

var_names = strvcat(xvars);
nyrs = endYear - begYear + 1;

%%
%==================================================================
%                       Estimate SDM Model
%==================================================================

[n,k] = size(x);
% k is all explanatory variables, including W-x

resulto = ols(yvec,x);
% vnames = strvcat('rnd',var_names);
% prt(resulto,vnames);


% 						Refer to Jim's textbook (equations 3.6 - 3.8) for elements
% 						of the algorithm. Specifically, instead of using
% 						 		(I-rho*W) y = XB + e
% 						it is easier to break the estimation into:
% 								y  = XBo + eo		and
%								Wy = XBd + ed
% 						so as not to have rho in the estimation.



rmin = -1;     % use -1,1 rho interval as default
rmax = 1;
sige = 1.0;
nu = 0;
d0 = 0;
c = zeros(k,1);   % diffuse prior for beta
T = eye(k)*1e+12;


% compute this stuff once to save time
TI = T\eye(k);
TIc = TI*c;


xpx = x'*x;
xpy = x'*yvec;

xpWy = x'*Wy;

bo = xpx\xpy;

bd = xpx\xpWy;

eo = yvec - x*bo;
ed = Wy - x*bd;

% ------------------------------------------- %
% Pre-calculate log-dets over a grid of rho values
rgrid = -0.99:0.01:0.99;
ngrid = size(rgrid,2);
detval = zeros(ngrid, 2);
time0 = clock;


% Start the parallel pool if not already running
if isempty(gcp('nocreate'))
    parpool('local', numWorkers); % Replace numWorkers with desired cores
end

hwait = waitbar(0,'calculating log-determinant ...');


parfor i = 1:ngrid
    rho = rgrid(i);

    detm = 0;
    for j = 1:nyrs

        % Calculate |I - rho*W|
        B = speye(nvec(j)) - rho * W{j};

        % Compute incomplete LU factorization
        [~, U] = ilu(B, struct('type','nofill')); % Obtain U factor

        % Extract diagonal elements of U
        du = diag(U);

        % Accumulate the sum of log-determinants
        detm = detm + sum(log(abs(du))); % Use abs to avoid complex logs
    end
    waitbar(i/ngrid);

    detval(i, :) = [rho, detm];
end
close(hwait);

% Optionally, display progress
% completed = sum(progress);
% fprintf('Completed %d out of %d iterations.\n', completed, ngrid);

% Close the parallel pool if desired
% delete(gcp('nocreate'));
% ------------------------------------------- %
%DONE UP TO HERE

%% Interpolate a finer grid
fgrid = -0.99:0.0001:0.99;
lndet = interp1(detval(:,1),detval(:,2),fgrid','spline');
detval = [fgrid' lndet];

det_time = etime(clock,time0);

ndraw = 600;

% storage for draws
bsave = zeros(ndraw,k);
psave = zeros(ndraw,1);
ssave = zeros(ndraw,1);
acc_rate = zeros(ndraw,1);
cc_save = zeros(ndraw,1);
cc = 0.1;
acc = 0;
gg=1; % for Griddy Gibbs sampling of rho
%     gg=0; % for Metropolis-Hastings sampling of rho
iter = 1;

rho = 0.1; % starting value

time1 = clock;

while (iter <= ndraw);
          % update beta
          AI = (xpx + sige*TI)\eye(k);
          ys = yvec - rho*Wy;
          b = x'*ys + sige*TIc;
          b0 = (xpx + sige*TI)\b;
          bhat = norm_rnd(sige*AI) + b0;
          xb = x*bhat;

          bsave(iter,:) = bhat';

          % update sige
          nu1 = n + 2*nu;
          e = (ys - xb);
          d1 = 2*d0 + e'*e;
          chi = chis_rnd(1,nu1);
          sige = d1/chi;

          ssave(iter,1) = sige;

		  % update rho using Griddy Gibbs
          if (gg == 1)

			  % Using the equations y  = XBo + eo and Wy = XBd + ed
			  %	get the betas and error terms e.
              b0 = (xpx + sige*TI)\(xpy + sige*TIc);
              bd = (xpx + sige*TI)\(xpWy + sige*TIc);
              e0 = yvec - x*b0;
              ed = Wy - x*bd;
              epe0 = e0'*e0;
              eped = ed'*ed;
              epe0d = ed'*e0;

              nmk = n/2;
              nrho = length(detval(:,1));
              iota = ones(nrho,1);

			  % For each potential value of rho, calculate the log likelihood.
              z = epe0*iota - 2*detval(:,1)*epe0d + detval(:,1).*detval(:,1)*eped;
              z = -nmk*log(z);
              den =  detval(:,2) + z;

              nn = length(den);
              yy = detval(:,1);

			  % Calculation of log likelihood missing a scaling value,
			  % and need an adjustment.
              adj = max(den);
              den = den - adj;
              xx = exp(den);

              % Use the trapezoid rule to integrate the area under the pdf.
              isum = sum((yy(2:nn,1) + yy(1:nn-1,1)).*(xx(2:nn,1) - xx(1:nn-1,1))/2);
              z = abs(xx/isum);

			  % Store the cdf.
              den = cumsum(z);

			  % Using a uniform distribution, randomly select a value of
			  % rho using the cdf (by inverting it).
              rnd = unif_rnd(1,0,1)*sum(z);
              ind = find(den <= rnd);
              idraw = max(ind);
              if (idraw > 0 && idraw < nrho)
                  rho = detval(idraw,1);
              end;

              psave(iter,1) = rho;

% GG ends here
          elseif (gg == 0)
              % update rho using Metropolis-Hastings
              % lookup log-det value for this value of rho
              gsize = detval(2,1) - detval(1,1);

              i1 = find(detval(:,1) <= rho + gsize);
              i2 = find(detval(:,1) <= rho - gsize);
              i1 = max(i1);
              i2 = max(i2);
              index = round((i1+i2)/2);

              if isempty(index)
                  index = 1;
              end;

              detm = detval(index,2);

              e = eo - rho*ed;
              epe = e'*e;
              epe = epe/(2*sige);
              rhox =  detm - epe - (n/2)*log(sige) - (n/2)*log(2*pi);

              rho2 = rho + cc*randn(1,1); % proposal for rho

              accept = 0;
              while accept == 0
                  if ((rho2 > rmin) && (rho2 < rmax));
                      accept = 1;
                  else
                      rho2 = rho + cc*randn(1,1);
                  end;
              end;

              % lookup log-det value for this value of rho
              gsize = detval(2,1) - detval(1,1);
              i1 = find(detval(:,1) <= rho2 + gsize);
              i2 = find(detval(:,1) <= rho2 - gsize);
              i1 = max(i1);
              i2 = max(i2);
              index = round((i1+i2)/2);
              if isempty(index)
                  index = 1;
              end;

              detm = detval(index,2);

              e = eo - rho2*ed;
              epe = e'*e;
              epe = epe/(2*sige);
              rhoy =  detm - epe - (n/2)*log(sige) - (n/2)*log(2*pi);

              % Metropolis-Hastings accept/reject step
              ru = unif_rnd(1,0,1);
              if ((rhoy - rhox) > exp(1)),
                  p = 1;
              else
                  ratio = exp(rhoy-rhox);
                  p = min(1,ratio);
              end;

              if (ru < p)
                  rho = rho2;
                  acc = acc + 1;
              end;

              psave(iter,1) = rho;

              acc_rate(iter,1) = acc/iter;

              % update cc based on std of rho draws
              if acc_rate(iter,1) < 0.4
                  cc = cc/1.1;
              end;
              if acc_rate(iter,1) > 0.6
                  cc = cc*1.1;
              end;

              cc_save(iter,1) = cc;

              if cc > 1
                  cc = 0.1;
              elseif cc < 1e-04;
                  cc = 0.1;
              end;
              %
          end;

		  % Plot the parameters from each draw.
          tt=1:iter;
          subplot(2,2,1),
          plot(tt,bsave(1:iter,1));
          xlabel('beta1');
          subplot(2,2,2),
          plot(tt,psave(1:iter,1));
          xlabel('rho');
          subplot(2,2,3),
          plot(tt,ssave(1:iter,:));
          xlabel('sige');
          drawnow;


          iter = iter+1;

end;

mcmc_time = etime(clock,time1);

% fprintf(1,'log-determinant time = %16.8f \n',det_time);
% fprintf(1,'MCMC sampling time = %16.8f \n',mcmc_time);
% fprintf(1,'ndraws = %10d \n',ndraw);

% discard the first observations
nomit = 500;

pdraws = psave(nomit+1:ndraw,1);

bdraws = bsave(nomit+1:ndraw,:);

sdraws = ssave(nomit+1:ndraw,1);

ng = length(sdraws);


%==================================================================
%                       Report the results
%==================================================================


% tt=1:ng;
% subplot(2,1,1),
% plot(tt,pdraws,'.');
% xlabel('rho draws');
% subplot(2,1,2),
% plot(tt,sdraws,'.');
% xlabel('sige draws');
% pause;
%
% subplot(2,2,1),
% plot(tt,bdraws(:,1),'.');
% xlabel(xnames(1,:));
% subplot(2,2,2),
% plot(tt,bdraws(:,2),'.');
% xlabel(xnames(2,:));
% subplot(2,2,3),
% plot(tt,bdraws(:,3),'.');
% xlabel(xnames(3,:));
% subplot(2,2,4),
% plot(tt,bdraws(:,4),'.');
%  xlabel(xnames(4,:));
% pause;

% subplot(1,1,1),
% plot(acc_rate,'.');
% pause;
% %
subplot(1,1,1),
[h1,f1,y1] = pltdens(pdraws);
plot(y1,f1,'.');
% xlabel('posterior density of rho');
%
pmean = mean(pdraws);
bmean = mean(bdraws);
smean = mean(sdraws);
% compute R-squared etc.
e = eo - pmean*ed;
epe = e'*e;
sigma = epe/(n-k);
ym = yvec - mean(yvec);
rsqr1 = epe;
rsqr2 = ym'*ym;
rsqr = 1- rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(n-k);
rsqr2 = rsqr2/(n-1);
rbar = 1 - (rsqr1/rsqr2); % rbar-squared
%%
fid = 1;
% fprintf(fid,'Spatial Durbin Model Estimates \n');
% % fprintf(fid,'Dependent Variable = %16s \n',vnames(1,:));
% fprintf(fid,'R-squared          = %9.4f \n',rsqr);
% fprintf(fid,'Rbar-squared       = %9.4f \n',rbar);
% fprintf(fid,'sigma^2 =epe/n-k   = %9.6f \n',sigma);
% fprintf(fid,'sige^2 posterior   = %9.6f \n',smean);
% fprintf(fid,'Nobs, Nvars        = %6d,%6d \n',n,k);
% fprintf(fid,'ndraws,nomit       = %6d,%6d \n',ndraw,nomit);
% fprintf(fid,'min and max rho    = %9.4f,%9.4f \n',rmin,rmax);
% % print timing information
% fprintf(fid,'MCMC time in secs  = %9.4f \n',mcmc_time);
% fprintf(fid,'time for lndet     = %9.4f \n',det_time);

results.rsqr		= rsqr		;
results.rbar		= rbar		;
results.sigma		= sigma		;
results.smean		= smean		;
results.n			= n			;
results.k			= k			;
results.ndraw		= ndraw		;
results.nomit		= nomit		;
results.rmin		= rmin		;
results.rmax		= rmax		;
results.mcmc_time	= mcmc_time	;
results.det_time	= det_time	;


tmp = pdraws;
bounds = cr_interval(tmp,0.95);
upper = bounds(:,1)';
lower = bounds(:,2)';


pout = [lower pmean upper];

bounds = cr_interval(tmp,0.99);
upper = bounds(:,1)';
lower = bounds(:,2)';

pout = [lower pout upper];

cnames = strvcat('lower01','lower05','mean','upper95','upper99');
rnames = strvcat('parameters','rho');


tmp = bdraws;
bounds = cr_interval(tmp,0.95);
upper = bounds(:,1);
lower = bounds(:,2);

bout =[lower bmean' upper];

bounds = cr_interval(tmp,0.99);
upper = bounds(:,1);
lower = bounds(:,2);

bout =[lower bout upper];

rnames = strvcat(rnames,var_names);

tmp = sdraws;
bounds = cr_interval(tmp,0.95);
upper = bounds(:,1);
lower = bounds(:,2);

sout = [lower smean upper];

bounds = cr_interval(tmp,0.99);
upper = bounds(:,1);
lower = bounds(:,2);

sout = [lower sout upper];

rnames = strvcat(rnames,'sige');

out = [pout
       bout
       sout];

in.rnames = rnames;
in.cnames = cnames;
in.fmt = '%16.6f';

% mprint(out,in);

results.rho_out = out;
results.rho_in = in;


bstring = 'Coefficient';
tstring = 't_stat';
pstring = 't_prob';
estring = 'std_err';
lstring = 'lower_01';
ustring = 'upper_99';
cnames = strvcat(bstring,tstring,pstring,estring,lstring,ustring);

cmean = mean([pdraws, bdraws, sdraws])';
smean = std([pdraws, bdraws, sdraws])';
bounds = cr_interval([pdraws, bdraws, sdraws],0.99);
ubounds = bounds(:,1);
lbounds = bounds(:,2);

results.b_out = [cmean cmean./smean tdis_prb(cmean./smean,sum(nvec)) smean lbounds ubounds];
results.b_in.rnames = rnames;
results.b_in.cnames = cnames;

%% calculate partial derivative effects estimates

% traces for the x-impacts calculations
% loop over all W-matrices

uiter=50;
maxorderu=100;
tracew=zeros(maxorderu,1);
traces = zeros(maxorderu,1);

for j=1:nyrs;

    nobs = nvec(j,1);
    rv=randn(nobs,uiter);

    wjjju=rv;
    for jjj=1:maxorderu
        wjjju=W{j}*wjjju;
        tracew(jjj)=mean(mean(rv.*wjjju));

    end

    % Directly calculate the (average) trace of W^2
    tracew(2)= sum(sum(W{j}'.*W{j}))/nvec(j,1);
    traces= traces + tracew;
    traces(1)=0;
    % traces(2)= sum(sum(W{j}'.*W{j}))/nvec(j,1);

end;

trs=[1;1/nyrs*traces];
ntrs=length(trs);
trbig=trs';


%% =============================================
% Below is a correction for SDM effects
trbig2 = [trbig(1,2:end) trbig(1,end)];
trmat = [trbig
        trbig2];
% =============================================


ree = 0:1:ntrs-1;

% we only want estimates for the non-fixed effects variables

rmat     = zeros(1,ntrs);
total    = zeros(ndraw-nomit,p,ntrs);
direct   = zeros(ndraw-nomit,p,ntrs);
indirect = zeros(ndraw-nomit,p,ntrs);

for iter=1:ndraw-nomit;

    %     rmat = zeros(1,ntrs);
    rmat = pdraws(iter,1).^ree;
    for j=1:p

        % beta = [bdraws(iter,j) bdraws(iter,j+p)];
        total(iter,j,:) = bdraws(iter,j)*rmat;

        direct(iter,j,:) = (bdraws(iter,j)*trbig).*rmat;
        indirect(iter,j,:) = total(iter,j,:) - direct(iter,j,:);
    end

end;

% cumulate over 101 orders
total_out = zeros(p,6);
total_save = zeros(ndraw-nomit,p);
for i=1:p;
    tmp = squeeze(total(:,i,:)); % an ndraw by 1 by ntraces matrix
    total_mean = mean(tmp);
    total_std = std(tmp);
    % Bayesian 0.99 credible intervals
    % for the cumulative total effects
    total_sum = (sum(tmp'))'; % an ndraw by 1 vector
    cum_mean = cumsum(mean(tmp));
    cum_std = cumsum(std(tmp));
    total_save(:,i) = total_sum;
    bounds = cr_interval(total_sum,0.99);
    cmean = mean(total_sum);
    smean = std(total_sum);
    ubounds = bounds(1,1);
    lbounds = bounds(1,2);
    total_out(i,:) = [cmean cmean./smean tdis_prb(cmean./smean,nobs) smean lbounds ubounds];
end;

% now do indirect effects
indirect_out = zeros(p,6);
indirect_save = zeros(ndraw-nomit,p);
for i=1:p;
    tmp = squeeze(indirect(:,i,:)); % an ndraw by 1 by ntraces matrix
    indirect_mean = mean(tmp);
    indirect_std = std(tmp);
    % Bayesian 0.95 credible intervals
    % for the cumulative indirect effects
    indirect_sum = (sum(tmp'))'; % an ndraw by 1 vector
    cum_mean = cumsum(mean(tmp));
    cum_std = cumsum(std(tmp));
    indirect_save(:,i) = indirect_sum;
    bounds = cr_interval(indirect_sum,0.99);
    cmean = mean(indirect_sum);
    smean = std(indirect_sum);
    ubounds = bounds(1,1);
    lbounds = bounds(1,2);
    indirect_out(i,:) = [cmean cmean./smean tdis_prb(cmean./smean,nobs) smean lbounds ubounds  ];
end;



% now do direct effects
direct_out = zeros(p,6);
direct_save = zeros(ndraw-nomit,p);
for i=1:p;
tmp = squeeze(direct(:,i,:)); % an ndraw by 1 by ntraces matrix
direct_mean = mean(tmp);
direct_std = std(tmp);
% Bayesian 0.95 credible intervals
% for the cumulative direct effects
direct_sum = (sum(tmp'))'; % an ndraw by 1 vector
cum_mean = cumsum(mean(tmp));
cum_std = cumsum(std(tmp));
direct_save(:,i) = direct_sum;
bounds = cr_interval(direct_sum,0.99);
cmean = mean(direct_sum);
smean = std(direct_sum);
ubounds = bounds(1,1);
lbounds = bounds(1,2);
direct_out(i,:) = [cmean cmean./smean tdis_prb(cmean./smean,nobs) smean lbounds ubounds  ];
end;


% now print x-effects estimates

bstring = 'Coefficient';
tstring = 't_stat';
pstring = 't_prob';
estring = 'std_err';
lstring = 'lower_01';
ustring = 'upper_99';
cnames = strvcat(bstring,tstring,pstring,estring,lstring,ustring);
ini.cnames = cnames;
ini.width = 2000;

% print effects estimates
vnameso = strvcat(xvars);
ini.rnames = strvcat('Direct  ',vnameso);
ini.fmt = '%16.6f';

% set up print out matrix
printout = direct_out;
% mprint(printout,ini);

results.direct_out = printout;
results.direct_in = ini;

printout = indirect_out;
ini.rnames = strvcat('Indirect',vnameso);
% mprint(printout,ini);

results.indirect_out = printout;
results.indirect_in = ini;

printout = total_out;
ini.rnames = strvcat('Total   ',vnameso);
% mprint(printout,ini);


results.total_out = printout;
results.total_in = ini;

results.eps = [ID (eo - pmean * ed)];

results.total    = total   ;
results.direct   = direct  ;
results.indirect = indirect;
results.bdraw = bdraws ;
results.pdraw = pdraws ;
results.W = W;
results.IDS = IDS;

